## [1] "Run Completed at 2015-12-05 13:33:19"

1 Aim

Which birds use which plants? What are the mechanisms that determine niche overlap in this assemblage.

We will fit several difference types of models and ask how they inform our understanding. For each model, we try to calculate a measure of discrepancy and assess model fit.

1.1 Data

  • Using automated cameras and transects, we have been recording flower visitation rates by each hummingbird species for two years.
  • We have bill and corolla lengths for each species of bird and flower.
  • We have been collecting phenology data on available resources

Birds use flowers that have similar length of bill shape and we expect the importance of trait matching to decrease as available resources increase.

2 Read in Data

## character(0)
##  [1] "Missing Trait Information:  sp"                       
##  [2] "Missing Trait Information: Burmeistera cylindrocarpa" 
##  [3] "Missing Trait Information: Burmeistera multiflora"    
##  [4] "Missing Trait Information: Calathea ischnosiphonoides"
##  [5] "Missing Trait Information: Canna jaegeriana"          
##  [6] "Missing Trait Information: Cavendishia tarapotana"    
##  [7] "Missing Trait Information: Drymonia brochidodroma"    
##  [8] "Missing Trait Information: Drymonia collegarum"       
##  [9] "Missing Trait Information: Gasteranthus pansamalanus" 
## [10] "Missing Trait Information: Gasternathus lateralis"    
## [11] "Missing Trait Information: Guzmania lehmanniana"      
## [12] "Missing Trait Information: Guzmania squarrosa"        
## [13] "Missing Trait Information: Justicia secunda"          
## [14] "Missing Trait Information: Macleania recumbens"       
## [15] "Missing Trait Information: Palicourea sodiroi"        
## [16] "Missing Trait Information: Pitcairnia nigra"          
## [17] "Missing Trait Information: Psamisia ???"              
## [18] "Missing Trait Information: Psammisia ecuadorensis"    
## [19] "Missing Trait Information: Stromantheus stromanthe"   
## [20] "Missing Trait Information: Zingibiraceae sp."

2.1 Match Species to Morphology

2.1.1 Elevation ranges

Create a binary variable whether each observation was in a low elevation or high elevation transect. We have some species that just occur at the top of the gradient, and are not present in the sampling window of flowers at the low elevation.

Accounting for non-availability. We have to figure out which plants were sampled in which periods, and if it was sampled, the non-detection are 0 if it wasn’t the non-detection are NA. then remove all the Na’s.

##                 Hummingbird  Low        m High Index
## 1   Green-crowned Woodnymph 1359 1376.645 1390     1
## 2 Rufous-tailed Hummingbird 1368 1476.333 1380     1
## 3            Andean Emerald 1378 1386.261 1380     1
## 4      White-necked Jacobin 1368 1397.500 1422     1
## 5    Stripe-throated Hermit 1368 1440.887 1450     1
## 6    White-whiskered Hermit 1340 1412.391 1450     1

3 Available Resources

What is the temporal distribution of resources?

Appears to be a log normal distribution with the existance of rare, but intense, pulses.

3.0.1 Summarize Observations

## Source: local data frame [1,301 x 7]
## Groups: Bird, Plant, ID [1204]
## 
##     Bird Plant     ID      DateP  Yobs  Elev Transect_R
##    (int) (int)  (chr)     (fctr) (int) (dbl)     (fctr)
## 1     15    85  FH709 2014-02-27    21  1990         NA
## 2     22    93  FH616 2014-05-07    17  2450         NA
## 3     17    85  FH414 2014-04-22    13  1990         NA
## 4     24   101  FL083 2013-07-29    13  2350         NA
## 5      6    98  FL061 2013-07-10    12  1325         NA
## 6     15    18  FH303 2013-11-19    12  1940         NA
## 7     15    85 FH1213 2015-02-11    11  1990         NA
## 8     24   101  FL083 2013-07-28    11  2350         NA
## 9      5     3  FL057 2013-07-05    10  1390         NA
## 10     8    92  NF084 2014-05-24    10  1513         NA
## ..   ...   ...    ...        ...   ...   ...        ...

What elevation transect is each observation in? The camera data need to be inferred from the GPS point.

4 Choosing a response distribution

Before we account for elevation and availability, what is the distribution of intensity of interactions?

If we tend to think of the intensity of events as poisson, our data is extremely ones-inflated. The first moment of the poisson does not well capture the data.

The best fitting moment is in opaque red. The observed data is in black.

What about negative binomial?

The best fitting moment is in opaque red. The observed data is in black.

What about log normal?

The best fitting moment is in opaque red. The observed data is in black.

5 Absences - accounting for non-detection

We have more information than just the presences, given species elevation ranges, we have absences as well. Absences are birds that occur at the elevation of the plant sample, but were not recorded feeding on the flower.

Remerge data with predictors

6 Resource availability at each data point

Three measures of resources.

Give binary versions of these resources as well - what matters is if its high or low for that resource measure. This helps account for the enormously long tail of resource pulses and minimizes the noise associated with such a widely ranging variable. The points on the huge blooms are going to have large leverage.

View quantile thresholds (0.5,.75,.9) as vertical lines.

6.0.1 How many counts for each quantile class?

Do we have enough data?

Relationship among resource measures

How many data points for each species for each resource level?

7 Exploratory Analysis

Let’s look at the same figures as above, with the addition of informed 0’s for non-detections.

7.1 Mean Intensity as a function of traits

7.2 Mean Intensity as a function of resources

Boxplots of trait-matching

7.3 Mixed Effect models

For initial inference, let’s look at some reasonable mixed effect models. There is no zero inflation here, so we treat these data with some caution.

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Yobs ~ Traitmatch + (1 | Bird)
##    Data: obsf
## 
##      AIC      BIC   logLik deviance df.resid 
##  13991.3  14014.0  -6992.6  13985.3    14624 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
##  -0.943  -0.441  -0.286  -0.164 101.435 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Bird   (Intercept) 1.064    1.031   
## Number of obs: 14627, groups:  Bird, 14
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.41516    0.27812  -5.088 3.61e-07 ***
## Traitmatch  -0.71907    0.02997 -23.993  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## Traitmatch -0.077
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Yobs ~ Traitmatch * BAll_Flowers + (1 | Bird)
##    Data: obsf
## 
##      AIC      BIC   logLik deviance df.resid 
##  13980.4  14018.4  -6985.2  13970.4    14622 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
##  -0.937  -0.440  -0.281  -0.163 103.570 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Bird   (Intercept) 1.06     1.03    
## Number of obs: 14627, groups:  Bird, 14
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -1.42350    0.36809  -3.867  0.00011 ***
## Traitmatch              -0.31884    0.15155  -2.104  0.03539 *  
## BAll_Flowers             0.01151    0.24729   0.047  0.96289    
## Traitmatch:BAll_Flowers -0.41433    0.15403  -2.690  0.00715 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmtc BAll_F
## Traitmatch  -0.535              
## BAll_Flowrs -0.659  0.799       
## Trtmtc:BA_F  0.527 -0.980 -0.804
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Yobs ~ Traitmatch * BUsed_Flowers + (1 | Bird)
##    Data: obsf
## 
##      AIC      BIC   logLik deviance df.resid 
##  13993.8  14031.7  -6991.9  13983.8    14622 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
##  -0.970  -0.442  -0.280  -0.163 100.585 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Bird   (Intercept) 1.068    1.033   
## Number of obs: 14627, groups:  Bird, 14
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -1.43827    0.27877  -5.159 2.48e-07 ***
## Traitmatch               -0.70505    0.03372 -20.909  < 2e-16 ***
## BUsed_Flowers             0.08064    0.06542   1.233    0.218    
## Traitmatch:BUsed_Flowers -0.05101    0.06144  -0.830    0.406    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmtc BUsd_F
## Traitmatch  -0.091              
## BUsed_Flwrs -0.068  0.346       
## Trtmtc:BU_F  0.045 -0.461 -0.685
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Yobs ~ Traitmatch * BFlowerA + (1 | Bird)
##    Data: obsf
## 
##      AIC      BIC   logLik deviance df.resid 
##  13920.7  13958.7  -6955.4  13910.7    14622 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -1.019 -0.459 -0.286 -0.166 91.558 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Bird   (Intercept) 1.128    1.062   
## Number of obs: 14627, groups:  Bird, 14
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.39891    0.28566  -4.897 9.73e-07 ***
## Traitmatch          -0.65461    0.03112 -21.033  < 2e-16 ***
## BFlowerA            -0.01816    0.06687  -0.272    0.786    
## Traitmatch:BFlowerA -0.42822    0.07371  -5.810 6.25e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmtc BFlwrA
## Traitmatch  -0.085              
## BFlowerA    -0.059  0.321       
## Trtmtch:BFA  0.034 -0.345 -0.661
## Data: obsf
## Models:
## tr_interactionA: Yobs ~ Traitmatch * BAll_Flowers + (1 | Bird)
## tr_interactionU: Yobs ~ Traitmatch * BUsed_Flowers + (1 | Bird)
## tr_interactionF: Yobs ~ Traitmatch * BFlowerA + (1 | Bird)
##                 Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## tr_interactionA  5 13980 14018 -6985.2    13970                         
## tr_interactionU  5 13994 14032 -6991.9    13984  0.000      0          1
## tr_interactionF  5 13921 13959 -6955.4    13911 73.055      0     <2e-16
##                    
## tr_interactionA    
## tr_interactionU    
## tr_interactionF ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

What does the mixed effect model look like?

7.3.1 Log tranformed presence only records

## Linear mixed model fit by REML ['lmerMod']
## Formula: log(Yobs) ~ Traitmatch * FlowerA + (1 | Bird)
##    Data: obsf[obsf$Yobs > 0, ]
## 
## REML criterion at convergence: 2338.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0024 -0.6775 -0.5439  0.5318  4.6726 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bird     (Intercept) 0.006485 0.08053 
##  Residual             0.339361 0.58255 
## Number of obs: 1296, groups:  Bird, 14
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)         3.968e-01  3.486e-02  11.383
## Traitmatch         -2.130e-02  2.212e-02  -0.963
## FlowerA             1.764e-05  1.247e-05   1.415
## Traitmatch:FlowerA -2.267e-05  1.479e-05  -1.534
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmtc FlowrA
## Traitmatch  -0.528              
## FlowerA     -0.274  0.261       
## Trtmtch:FlA  0.131 -0.246 -0.613
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

8 Hierarchical Bayesian Model

8.1 Model Equation

\[ Yobs_{i,j,k} \sim Binom(Detect_i,N_{i,j,k}) \]

\[ N_{i,j,k} \sim Pois(\lambda_{i,j,k}) \]

\[log(\lambda_{i,j,k})<-\alpha_i + \beta_i * abs(Bill_i - Corolla_j) * \beta_2 * Resources_k + \beta_3 * Resources_k * abs(Bill_i - Corolla_j) \]

Priors

\[ Detect \sim N(dprior,tau_detect) \] \[\alpha_i \sim N(intercept,\tau_{\alpha})\] \[\beta_{1,i} \sim N(\gamma_1i,\tau_{\beta_1})\] \[\beta_{2,i} \sim N(\gamma_2i,\tau_{\beta_2})\] \[\beta_{3,i} \sim N(\gamma_3i,\tau_{\beta_3})\]

Hyperpriors

Group Level Means \[dprior \sim U(0,0.5)\] \[\gamma_{1,i} \sim N(0,0.0001)\] \[\gamma_{2,i} \sim N(0,0.0001)\] \[\gamma_{3,i} \sim N(0,0.0001)\] \[ intercept \sim N(0,0.0001)\]

Group Level Variance

\[\tau_{\alpha} \sim Gamma(0.0001,0.0001)\] \[\tau_detect \sim Gamma(0.0001,0.0001)\] \[\tau_\beta1 \sim Gamma(0.0001,0.0001)\] \[\tau_\beta2 \sim Gamma(0.0001,0.0001)\] \[\tau_\beta3 \sim Gamma(0.0001,0.0001)\]

Derived quantities

\[\sigma_{int} = \sqrt[2]{\frac{1}{\tau_\alpha}}\] \[\sigma_detect = \sqrt[2]{\frac{1}{\tau_detect}\] \[\sigma_{slope1} = \sqrt[2]{\frac{1}{\tau_\beta1}}\] \[\sigma_{slope2} = \sqrt[2]{\frac{1}{\tau_\beta2}}\] \[\sigma_{slope3} = \sqrt[2]{\frac{1}{\tau_\beta3}}\]

## 
## sink("Bayesian/NmixturePoissonRagged.jags")
## 
## cat("
##     model {
##     for (x in 1:Nobs){
##     
##     # Covariates for true state   
##     log(lambda[Bird[x],Plant[x],Time[x]]) <- alpha[Bird[x]] + beta1[Bird[x]] * traitmatch[x] + beta2[Bird[x]] * resources[x] + beta3[Bird[x]] * resources[x] * traitmatch[x]      
##     
##     #True State
##     N[x] ~ dpois(lambda[Bird[x],Plant[x],Time[x]] )    
##     
##     #Observation Process
##     Yobs[x] ~ dbin(detect[Bird[x]],N[x])    
## 
##     #Assess Model Fit
##       
##     #Fit discrepancy statistics
##     eval[x]<-detect[Bird[x]]*N[x]
##     E[x]<-pow((Yobs[x]-eval[x]),2)/(eval[x]+0.5)
## 
##     ynew[x]~dbin(detect[Bird[x]],N[x])
##     E.new[x]<-pow((ynew[x]-eval[x]),2)/(eval[x]+0.5)
## 
##     }
##     
##     for (i in 1:Birds){
##     detect[i] ~ dnorm(dprior,tau_detect)
##     alpha[i] ~ dnorm(intercept,tau_alpha)
##     beta1[i] ~ dnorm(gamma1,tau_beta1)    
##     beta2[i] ~ dnorm(gamma2,tau_beta2)    
##     beta3[i] ~ dnorm(gamma3,tau_beta3)    
##     }
##   
## 
##     #Hyperpriors
##     #Slope grouping
##     gamma1~dnorm(0,0.0001)
##     gamma2~dnorm(0,0.0001)
##     gamma3~dnorm(0,0.0001)
##     
##     #Intercept grouping
##     intercept~dnorm(0,0.0001)
##   
##     #detect grouping
##     dprior~dunif(0,0.5)
##     
##     # Group intercept variance
##     tau_alpha ~ dgamma(0.0001,0.0001)
##     sigma_int<-pow(1/tau_alpha,0.5) #Derived Quantity
##     
##     #Group Slope Variance
##     tau_beta1 ~ dgamma(0.0001,0.0001)
##     tau_beta2 ~ dgamma(0.0001,0.0001)
##     tau_beta3 ~ dgamma(0.0001,0.0001)
##     
##     sigma_slope1<-pow(1/tau_beta1,0.5)
##     sigma_slope2<-pow(1/tau_beta2,0.5)
##     sigma_slope3<-pow(1/tau_beta3,0.5)
##    
##     #Group Detection variance
##     tau_detect ~ dgamma(0.0001,0.0001)
##     sigma_detect<-pow(1/tau_detect,0.5)
## 
##     #derived posterior check
##     fit<-sum(E[]) #Discrepancy for the observed data
##     fitnew<-sum(E.new[])
##     
##     }
##     ",fill=TRUE)
## 
## sink()
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 237064
## 
## Initializing model

Re-update if needed to extend

8.2 Assess Convergence

8.3 Posteriors

8.4 Predicted Relationship

Hold resource relationship at 0, what does trait-matching look like?

8.5 Full posterior prediction

8.6 Trait-matching at High and Low Resource Density

That looks pretty crazy, the posterior is very broad, what about holding at mean value for all coefficients

Much cleaner, more iterations needed?

8.7 Visualize interactions

8.8 Species Predictions

8.8.1 Species Predictions: High and Low Resolutions

8.9 Species Level Interaction

8.10 Interaction density functions

Let’s take a closer look at distribution of interaction effect posteriors values for each species.

8.11 Interaction and Bill Length

Do species with long bill lengths have positive interaction effects?

8.12 Detection Probability

Detection Table

##                Hummingbird mean lower upper
## 5  Fawn-breasted Brilliant  0.6   0.2   1.2
## 10    Speckled Hummingbird  1.1   0.7   1.9
## 9   Purple-bibbed Whitetip  1.6   0.5   3.1
## 8  Green-fronted Lancebill  1.8   1.1   2.8
## 6        Gorgeted Sunangel  2.3   1.6   3.7
## 4            Collared Inca  2.4   1.2   3.7
## 11  Stripe-throated Hermit  2.4   1.6   3.6
## 2               Brown Inca  2.5   1.4   3.8
## 3      Buff-tailed Coronet  2.8   1.5   5.2
## 13     Violet-tailed Sylph  2.9   2.0   3.6
## 12    Tawny-bellied Hermit  3.5   2.7   4.4
## 7  Green-crowned Woodnymph  3.7   2.3   5.4
## 14  White-whiskered Hermit  4.0   2.8   5.6
## 1       Booted Racket-tail  4.6   3.5   6.0

9 Model Evaluation

9.1 Correlation coefficient among posteriors

Hierarchical Posteriors

Species level posteriors

9.2 Posterior Check

Order predictions from worst to best.

##                Hummingbird             Iplant_Double     chisq Yobs
## 1     Speckled Hummingbird        Palicourea lineata 350.23827   17
## 2               Brown Inca Mezobromelia capituligera 260.72200   21
## 3        Gorgeted Sunangel         Psammisia sodiroi 172.01422   13
## 4               Brown Inca           Bomarea pardina 152.53815   12
## 5      Violet-tailed Sylph Mezobromelia capituligera 137.51879   13
## 6        Gorgeted Sunangel         Psammisia sodiroi 128.64470   11
## 7               Brown Inca        Heliconia impudica 105.08869   10
## 8            Collared Inca        Pitcairnia sodiroi  94.21488    9
## 9     Speckled Hummingbird         Macleania stricta  92.33434    8
## 10 Green-fronted Lancebill    Palicourea acetosoides  89.32503    8
##    Used_Flowers Traitmatch
## 1    -0.6396768      0.140
## 2    -0.5886063      0.156
## 3     1.9224631      0.408
## 4     0.1100659      1.988
## 5    -0.5349655      1.576
## 6     1.9224631      0.408
## 7     0.4031075      1.585
## 8    -0.6932043      1.576
## 9     0.1239308      0.536
## 10   -0.6532197      0.964